#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <omp.h>
#include <mpi.h>

#define TRUE    1;
#define FALSE   0;

/* Armazena os elementos para os cálculos do método Jacobi-Richardson*/
typedef struct {
    int J_ORDER;
    int J_ROW_TEST;
    float J_ERROR;
    int J_ITE_MAX;
    float **MA;
    float *MB;
} Jacobi;

/* Armazena o vetor com as soluções encontradas e o número de iterações */
typedef struct {
    float *S_SOLUTION;
    float *S_SOLUTION_AUX;
    int S_ITE;
} Solucao;

typedef struct {
	float valor;
	int id;
} valorID;

/* Faz a leitura dos valores necessários para os cálculos */
void inicializa_Jacobi(Jacobi *J, FILE *file){
    int i, j;

    fscanf(file, "%d", &(J->J_ORDER));
    fscanf(file, "%d", &(J->J_ROW_TEST));
    fscanf(file, "%f", &(J->J_ERROR));
    fscanf(file, "%d", &(J->J_ITE_MAX));

    J->MA = (float**) malloc(J->J_ORDER * sizeof(float*));
    for(i=0; i < J->J_ORDER; i++)
        J->MA[i] = (float*) malloc(J->J_ORDER * sizeof(float));

    J->MB = (float*) malloc(J->J_ORDER * sizeof(float));

    for(i=0; i < J->J_ORDER; i++){
        for(j=0; j < J->J_ORDER; j++){
            fscanf(file, "%f", &(J->MA[i][j]));
        }
    }
    #pragma omp for		
    for(i=0; i < J->J_ORDER; i++)
        fscanf(file, "%f", &(J->MB[i]));
}

/* Estritamente Diagonalmente Dominante: */
/* Verifica se a soma dos elementos de uma linha, retirando-se a diagonal, será menor que a diagonal */
/* Caso a soma seja igual o maior o método de Jacobi-Richardson irá divergir */
int converge(Jacobi *J){
    int i, j;
    int aux;

    for(i=0; i < J->J_ORDER; i++){
        aux = 0;

        for(j=0; j < J->J_ORDER; j++){
            if(i != j)
                aux += J->MA[i][j];
        }

        if(aux >= J->MA[i][i])
            return FALSE;
    }

    return TRUE;
}

/* Imprimi o vetor de soluções segundo o padrão requisitado */
void imprimi_Solucao(Jacobi J, Solucao S){
    int i;
    float result = 0.0f;
    #pragma omp parallel for reduction(+:result)	
    for(i=0; i < J.J_ORDER; i++)
        result += S.S_SOLUTION[i] * J.MA[J.J_ROW_TEST][i];
	
	printf("-------------------------------------------------------\n");
    printf("Iterations : %d\nRowTest: %d => [%f] ?= %f\n", S.S_ITE, J.J_ROW_TEST, result, J.MB[J.J_ROW_TEST]);
	printf("-------------------------------------------------------\n");
}

/* Método Jacobi-Richardson */
void Jacobi_Richardson(Jacobi J, Solucao *S, int id, int num_proc){
    	int i, j;
	int flag_fim = FALSE;
	int tam_calc;
	int namesize;
	char name[MPI_MAX_PROCESSOR_NAME];
	MPI_Status status;
	
	tam_calc = J.J_ORDER/(num_proc-1);
	MPI_Get_processor_name(name,&namesize);

	/* Calcula os resultados */
	if(id == 0){
		int ini;
		int fim;
		int w,z,y;
		float *valores_recebidos;
		float aux1 = 0.0f, aux2 = 0.0f, erro1 = 0.0f, erro2 = 0.0f;
		
	
		do{
			
			/* Envia de flag para não terminar os processos */
			MPI_Bcast(&flag_fim, 1, MPI_INT, 0, MPI_COMM_WORLD);

			/* Envio do vetor de soluções para os processos */
			MPI_Bcast(&(S->S_SOLUTION[0]), J.J_ORDER, MPI_FLOAT, 0, MPI_COMM_WORLD);
			
			for(i = 0 ; i < J.J_ORDER ; i++)
				S->S_SOLUTION_AUX[i] = S->S_SOLUTION[i];
			
			/* Recebimento dos novos valores */
			#pragma omp for private(j, w, y,ini, fim, valores_recebidos)		
		    	for(i=1; i < num_proc; i++){

				ini = (i-1)*tam_calc;
				if(i == num_proc-1)
					fim = J.J_ORDER-1;
				else
					fim = (i*tam_calc)-1;

				/*valores_recebidos = malloc((fim-ini+1)*sizeof(float));*/
				z = 0;
				
				MPI_Recv(valores_recebidos, (fim-ini+1), MPI_FLOAT,i, 0, MPI_COMM_WORLD, &status);
				
				for(w = ini ; w < fim+1 ; w++){
					S->S_SOLUTION[w] = valores_recebidos[z];
					z++;
				}
				erro1 = 0.0f;
		    		erro2 = 0.0f;
		    
				for(y=ini; y <= fim; y++){
		        		
					aux1 = S->S_SOLUTION[y] - S->S_SOLUTION_AUX[y];
					aux2 = S->S_SOLUTION[y];

					if(aux1 < 0)
					    aux1 *= -1;
					if(aux2 < 0)
					    aux2 *= -1;
					if(aux1 > erro1)
					    erro1 = aux1;
					if(aux2 > erro2)
					    erro2 = aux2;
					if((erro1/erro2) < J.J_ERROR)
						break;
		    		}

			}
		    S->S_ITE++;
	

		}while(((erro1/erro2) > J.J_ERROR) && (S->S_ITE < J.J_ITE_MAX));
		
		free(valores_recebidos);
		flag_fim = TRUE;

		/* Envio da flag para termino execução */
		MPI_Bcast(&flag_fim,1,MPI_INT,0,MPI_COMM_WORLD);
	}else{
		int ini;
		int fim;
		float result = 0.0f;
		float *valores;
			
		/*float aux1 = 0.0f, aux2 = 0.0f, erro1 = 0.0f, erro2 = 0.0f, erro_final = 1.0f;*/

		do{
			/* Recebe a mensagem do nó 0 */
			MPI_Bcast(&flag_fim, 1, MPI_INT, 0, MPI_COMM_WORLD);

			/* Caso a flag de fim recebida seja TRUE devido aos critérios de parada no nó 0,
			   finaliza o processo */			
			if(flag_fim == 1)
				break;

			/* Recebe o vetor de soluções */
			MPI_Bcast(&(S->S_SOLUTION[0]), J.J_ORDER, MPI_FLOAT, 0, MPI_COMM_WORLD);

			/*calculo das posiçoes de inicio e fim que devem ser calculadas por esse processo*/
			ini = (id-1)*tam_calc;
			if(id == num_proc-1)
				fim = J.J_ORDER-1;
			else
				fim = (id*tam_calc)-1;

			valores = malloc((fim-ini+1)*sizeof(float));
			

			#pragma omp parallel for private(j,result)
			for(i=ini; i <= fim; i++){
			        result = J.MB[i];

				for(j=0; j < J.J_ORDER; j++){
				    if(i != j)
				        result-= S->S_SOLUTION[j] * J.MA[i][j];
				}

				valores[i-ini] = result / J.MA[i][i];
			}

			MPI_Send(valores, (fim-ini+1), MPI_FLOAT, 0, 0, MPI_COMM_WORLD);

		    S->S_ITE++;

		}while(S->S_ITE < J.J_ITE_MAX);

		free(valores);
	}
	
}

int main(int argc, char **argv){
    Jacobi J;
    Solucao S;
    FILE *data;
	int id=0, p=0;

	/* Inicia os processos MPI */
	MPI_Init(&argc, &argv);
	MPI_Comm_rank(MPI_COMM_WORLD, &id);
	MPI_Comm_size(MPI_COMM_WORLD, &p);

    data = fopen("data-t3/matriz1000.txt","r");

    inicializa_Jacobi(&J, data);

    fclose(data);

	/* Aloca os vetores de solução */	
	S.S_ITE = 0;
    S.S_SOLUTION = (float*) malloc(J.J_ORDER * sizeof(float));
    S.S_SOLUTION_AUX = (float*) malloc(J.J_ORDER * sizeof(float));

	/* 
	 * Caso o sistema seja convergente, é aplicado o método Jacobi-Richardson
	 * caso contrário o sistema é divergente.
	*/    

	if(converge(&J)){
        Jacobi_Richardson(J, &S, id, p);
    }
    else{
		MPI_Abort(MPI_COMM_WORLD,-1);
		printf("Sistema linear divergente!");
		return EXIT_FAILURE;
	}
	
	if(id == 0)
		imprimi_Solucao(J, S);
	MPI_Finalize();

    return 0;
}

